Sensitivity and spectral control of network lasers

Recently, random lasing in complex networks has shown efficient lasing over more than 50 localised modes, promoted by multiple scattering over the underlying graph. If controlled, these network lasers can lead to fast-switching multifunctional light sources with synthesised spectrum. Here, we observe both in experiment and theory high sensitivity of the network laser spectrum to the spatial shape of the pump profile, with some modes for example increasing in intensity by 280% when switching off 7% of the pump beam. We solve the nonlinear equations within the steady state ab-initio laser theory (SALT) approximation over a graph and we show selective lasing of around 90% of the strongest intensity modes, effectively programming the spectrum of the lasing networks. In our experiments with polymer networks, this high sensitivity enables control of the lasing spectrum through non-uniform pump patterns. We propose the underlying complexity of the network modes as the key element behind efficient spectral control opening the way for the development of optical devices with wide impact for on-chip photonics for communication, sensing, and computation.


Supplementary Note 1: NetSALT: extending SALT theory to lasing networks
Here we provide a more detailed derivation of netSALT, the numerical model used in the simulations of lasing networks presented in this work. The accompanying code is available at https://github.com/arnaudon/netSALT.

A. Open quantum graphs
A quantum graph is a metric graph (i.e., a graph where each edge (ij) has an associated length l ij with an accompanying real length variable x ∈ [0, l ij ]) such that a function η(x), which is dened on each edge and thus dened on the entire graph, satises a dierential equation. In this case, we consider a quantum graph where the function η(x) satises the Helmholtz dierential equation where the complex-valued n ij ∈ C correspond to the index of refraction of the edge (ij).

This equation being linear, it has solutions of the form
where the complex-valued λ ± ij ∈ C represent the left-and right-propagating wave amplitudes. The continuity of η(x) at each node is ensured by considering the edge function η ij (x) evaluated at the nodes: η i := η ij (0) and η j := η ij (l ij ). ( Then the conservation of energy at each node i can be shown to be equivalent to where the sum is over the nodes adjacent to i, and the matrix L(k), dependent on the wavenumber, is a matrix acting on the node vector η with components η i dened in (3). We refer to [1,2] for more details.
The matrix L(k) can be expressed in terms of an extension of the graph incidence matrix, which allows the simplication of the calculations of various quantities (see Equation (27) and [3]). The condition (4) corresponds to an eigenvalue problem so, equivalently one can solve the corresponding scalar equation det (L(k)) = 0 , for discrete wavenumbers k µ . Numerically, we solve this equation using the smallest eigenvalue of L(k), which is ecient to compute with sparse matrices. Notice that the node representation contains a denominator term that diverges when kn ij l ij → nπ, with n = 1, 2, . . ..
This can cause numerical instabilities in rare cases, only encountered so far for graphs with several edges with the same length. We numerically x it by adding a small noise to the edge lengths (or node positions).
Edges with one open end (i.e., a node of degree 1) are considered to be outside of the cavity and admit no incoming wave. This can be simply written as a projection of the matrix L, where elements corresponding to the outgoing waves are projected out, thus allowing them to take any value (and not enforced to be vanishing from the right hand side of (4)). This condition makes the quantum graph open, or lossy, and any solution of (4) must have a complex wavenumber k µ . Henceforth, we make the distinction between inner edges, corresponding to the lasing cavity, and outer edges, corresponding to the open boundary of the cavity, also known as the last scattering surface in laser theory. For example, to describe the integration over the inner edges of the cavity, we use the shorthand notation For each passive mode k µ , the standard Q-factor is given by

B. The SALT equation on network edges
The SALT equation [4] describes the interaction of lasing modes under non-uniform pumping. In our formulation, the pump is naturally described by an amplitude D 0 and an edge indicator vector, denoted δ pump , with components δ pump,ij = 1 if the pump is illuminating edge (ij) and δ pump,ij = 0 otherwise. We will also use the notation δ pump (x) for the pump density, where δ pump (x) = δ pump,ij /l ij for x on edge (ij).
The lasing modes are the modes with Im(k µ ) = 0, and have the form where I µ is the modal intensity, and u µ is the mode prole, normalised as On edge (ij), the SALT equation [4] is a nonlinear extension of the Helmholtz equation (1) given by where γ µ = γ ⊥ kµ−ka+iγ ⊥ is a function that denes the gain spectrum with centre at k a and width γ ⊥ ; and Γ µ = −Im(γ µ ) is the Lorentzian gain curve.
We will not re-derive this equation here from several approximations of the Maxwell-Bloch equation, but refer to [4] for more details and only mention that key to the SALT model is the steady state assumption, or stationary inversion approximation, where the inversion population (denoted by D(x, t) in [5]) is taken to be constant in time. We refer to [57] for more detailed studies on the validity and generalisations of this approximation.

C. Finding threshold lasing modes
Before computing the modal intensities I µ , we need the threshold lasing modes, which are obtained as solutions of the linear equation This is an implicit equation for the lasing threshold D µ,th , the threshold wavenumber k µ,th , and the threshold lasing mode prole u µ,ij .
Solving this equation involves an iterative algorithm on the value of D 0 to reach the condition Im(k µ,th ) = 0, where the secular equation (6) is solved at each step. When D 0 is updated, we use the so-called Brownian ratchet algorithm [8] to search for the corresponding k µ (D 0 ). This algorithm proposes random moves in the complex plane of wavenumbers, and accepts only the ones that decrease the smallest eigenvalue of L(k), and stops the search when a certain threshold is reached. The size of the proposed moves is adjusted according to how far we expect the mode to have moved.
To speed up the search of threshold lasing modes, we rst estimate the location of a mode with a dierent D 0 by assuming that the mode proles do not change with the pump, i.e. u µ,ij (D 0 ) = η µ,ij . First, recall that η µ,ij are the passive modes, i.e., the solution of Multiplying (11) by η µ,ij and integrating over the cavity, we obtain where f µ is the pump overlapping factor of mode µ, dened as To estimate the pump strength at threshold, we use Im(k µ,th ) = 0 in (13) to get where we also used the fact that Real(γ µ ) is small for high Q modes.
To obtain an estimate of the complex wavenumber for an updated pump power D 0 = D 0 + δD 0 , we use, instead of (13), the equation where γ µ is now evaluated at k µ (D 0 ). This equation is obtained similarly to equation (13) by replacing the passive mode with a pumped mode.
To nd the threshold lasing modes, we linearly increase D 0 with small steps using (16) as a starting point for the Brownian ratchet algorithm to nd the next partially pumped mode until we reach Im(k) = 0, and then use a binary search (together with Brownian ratchet) to locate the position of the lasing threshold D µ,th .

D. Interacting modal intensities
Once the threshold lasing modes are found, we can estimate their modal intensities as a function of the pump power D 0 . To do this, we assume that the mode proles above threshold are the same as the mode proles at threshold, and the threshold wavenumbers k µ,th remain the same above threshold. With this approximation, which corresponds to the single pole approximation of [4], we can estimate the modal intensities of each mode, given a pump prole δ pump and a pump strength D 0 , as follows.
From (10), and using the normalisation (9), we follow [4] to arrive at the matrix equation where the sum is over lasing modes only, and the interaction matrix T has elements where the use of the real part is an approximation, since the integral in (18) has a small complex part in general.
Given D 0 , the modal intensities are then simply found as if the set of lasing modes (indexed as ν) are known. To nd the lasing modes, we follow again Ref. [4], and rst compute the interacting lasing thresholds D µ,int . For the rst lasing mode, the interaction threshold will be the lasing threshold D µ,th , but for the next lasing modes, interaction with the currently lasing modes will increase this value, until it reaches ∞, and no more modes can lase, leading to the phenomenon known as gain clamping.
To compute the sequence of lasing threshold modes, we proceed as follows. Let us assume that we have N lasing modes, and we seek to compute the interacting threshold of the next mode, indexed µ N +1 . At D 0 = D µ N +1 ,int , the mode µ N +1 will not lase, so I µ N +1 = 0, which, after some manipulation, gives an implicit equation for the interacting lasing threshold. Due to linearity, we can simply rearrange terms to get point, the denominator will become negative, corresponding to the gain clamping regime, where all other modes are suppressed by currently lasing modes, see [4] for more details on this phenomenon. Sometimes a lasing mode can stop lasing, due to a negative slope in (19), in which case this mode is removed from the list of lasing modes and will not contribute to this equation in the search of the next lasing mode.
The solution of this equation thus provides the so-called LL curves, with modal intensities of all the modes as a function of the pump power D 0 , given as piece-wise linear functions, as well as the possibility to approximate lasing spectra at a given pump power, if a lasing linewidth is added.

E. Pump optimisation in netSALT with linear programming
To numerically optimise the pump prole in netSALT, we do not evaluate the full modal intensities, as this would result in a costly and slow algorithm, due to the search of modal trajectories in the complex plane to reach threshold. Instead, we use the linear approximation of the lasing threshold (15), which depends on the pump overlapping factor f µ (δ pump ) given by (14). This overlapping factor can be written as a scalar product as it is given by a sum over the inside edges: where the component f µ,ij is the overlapping factor for a pump dened to be applied only to edge (ij). The optimal pump δ pump,µ is then the result of a minimisation problem of the form where a ν := f ν Q ν Γ ν and the hyper-parameter > 0 is a regulariser to avoid small pump proles. From this optimisation, we obtain a family of solutions with various coverages of the network surface area on the target mode prole.
The original optimisation (22) is combinatorial since the edge indicator vector δ is integervalued. To make the optimisation more amenable, we relax the problem to search for a real edge vector x, where each component (corresponding to an edge) is real and bounded To solve this relaxed problem, we rewrite it as the following linear program (LP): where m := max ν a T ν x and we have used the Charnes-Cooper transformation To solve the LP (24), we use the public python software PuLP available at https://github. com/coin-or/pulp.
From the continuous solution of the LP (23), we obtain a discrete pump vector, δ * µ , where δ * µ,ij = 1 if x ij > 0, followed by further thresholding of edges that have a small impact on the cost. Such edges are thought of as`noise' from the SALT approximation, which reduce the resulting modal suppression ratio.
The results of this optimisation on the Buon graph are illustrated in Supplementary The threshold lasing frequencies and the non-interacting lasing thresholds (not shown) also match the values reported in [4]. The code for this example is available in the github repository (https://github.com/arnaudon/netSALT).

B. Ring laser
In a ring with real index n and length L, the modes lie on the real axis with k m = 2mπ nL , where m is a positive integer. To model a ring laser with a nite Q-factor, we require complex refractive indexñ on the edges. Letñ = n + iκ. Then the modes are given by complex values: Loss is therefore dened via κ, or equivalently by the Q-factor (8). Supplementary Figure 13 shows the netSALT calculation of a uniformly pumped micro-ring laser with cavity length 10 µm and refractive index 1.5 + 0.005i. This example is also provided in the github repository (https://github.com/arnaudon/netSALT).

Supplementary Note 3: Additional calculations in netSALT
We collect here additional formulae of the netSALT model introduced above.
A. Matrix representation of L(k) The quantum graph equation (5) can be written in terms of matrices with analogues in classical graph theory. Given a network with N nodes and E edges, the matrix L(k) can be interpreted and rewritten as a quantum graph weighted Laplacian of the form where the matrix B(k) N ×2E is an extension of the standard incidence matrix in graph theory with elements which contains both edge directions, and the diagonal weight matrix W (k) 2E×2E is dened as W ij,ij = e 2ikl ij − 1 , and the same for W ji,ji . We refer to [3] for the details of the derivation of these equations.

B. Pump overlapping factor
The pump overlapping factor dened in (14) is explicitly given as where the block-diagonal matrix Z is formed by the blocks l ij e ikl ij l ij e ikl ij e 2ikl ij −1 2ik   .

C. Mode competition matrix
To simplify notation in this calculation, we consider k µ to be complex and incorporate the index of refraction and the pump term into γ, and we drop edge indices ij in the expressions (28) below.
The T matrix dened in (18) has elements and the integral can be written as

D. Calculation of the edge mean of |E| 2
The mode solution has the form (2): For brevity, we remove edge subscripts and take the modulus squared: If we integrate from x = 0 to l, we get This can be expressed in matrix form as which can be computed from the node solution. Here x stands for the complex conjugate of x.

E. Calculation of the Inverse Participation Ratio
The inverse participation ration (IPR) provides a measure for the mode spread over the graph, and is given by: where L tot is the total edge length. This formula can be evaluated analytically using the length of the un-pumped edges is 5% of the total graph length. Spectrum under uniform pumping (grey), also at D 0 = 0.02, is shown on the right panel for comparison. b-c The same calculations and plots are shown for two Buon graphs with low and high edge density, respectively. The Buon graphs clearly have higher spectral sensitivity compared to the single-loop graph. With increasing edge density, the total graph length increases corresponding to increase in the number of lasing modes, and larger spectral sensitivity. Density also aects the spatial extent of modes in the network, which in-turn has an eect on sensitivity as it changes the overlap between modes and the un-pumped edges. d Inverse participation ratio (IPR) distribution for the graphs in a-c shows that dense graphs (with more cycles) have larger fraction of localised modes and more number of competing lasing modes. The network lasers examined in the manuscript have edge densities between the low and dense Buon graphs shown in b and c.